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Confined qnantnm systems involving N identical interacting fermions are found in many areas of 
physics, including condensed matter, atomic, nuclear and chemical physics. In a previous series of 
papers, a manybody perturbation method that is applicable to both weakly and strongly-interacting 
systems of bosons has been set forth by the author and coworkers. A symmetry invariant pertur¬ 
bation theory was developed which uses group theory coupled with the dimension of space as the 
perturbation parameter to obtain an analytic correlated wave function through first order for a 
system under spherical confinement with a general two-body interaction. In the present paper, we 
extend this formalism to large systems of fermions, circumventing the numerical demands of ap¬ 
plying the Pauli principle by enforcing the Pauli principle on paper. The method does not scale in 
complexity with N and has minimal numerical cost. We apply the method to a unitary Fermi gas 


and compare to recent Monte Carlo values. 

PACS numbers: 03.65.Ge,03.75.Ss,31.15.xh,31.15.xp 

Introduction - Confined quantum systems of fermions 
are widespread across physics. They include, for ex¬ 
ample, atoms, atomic nuclei, neutron stars, quantum 
dots, and cold Fermi gases. These systems possess from 
a few tens to millions of particles and span an enor¬ 
mous range of interparticle interaction strength present¬ 
ing a challenge for Al-body methods when mean-field ap¬ 
proaches fail. In the last decade, ultracold Fermi gases 
have emerged as a testing ground for manybody methods 
due to their precise controllability in experiments. Using 
a magnetic field to tune the scattering length of atoms 
allows the exploration of the physics over many length 
scales, including the “unitary gas” defined by an infinite 
scattering length. Without a defining length scale other 
than the interparticle distance, the unitary gas does not 
yield to conventional perturbation treatments, typically 
requiring a full manybody treatment for an accurate de¬ 
scription. 

Fermi systems in the unitary regime are currently of 
great interest. This strongly interacting regime, stabi¬ 
lized by the Pauli exclusion principle, exists on the cusp 
of the BCS-BEC crossover and exhibits universal ther¬ 
modynamic behavior which has been verified in the lab¬ 
oratory to within a few percent [T| . Verifying this univer¬ 
sal behavior theoretically requires access to the partition 
function, i.e. the energy spectrum with degeneracies. Re¬ 
quired resources on a classical computer for an exact so¬ 
lution of even the ground state of the V-body problem 
scale exponentially with N. Traditional full configuration 
interaction can cope with at most some ten particles 0 
and methods such as coupled-cluster Q with a compu¬ 
tational time polynomial in N, 0{N'^), are extremely 
expensive. Quantum Monte Carlo (QMC) methods for 
fermions suffer from the so-called “sign” problem[3-@| 
resulting in an exponential growth in simulation times. 
Other methods that have been applied to strongly inter¬ 
acting Fermi systems include the method of correlated 


basis functions 0, [1], density functional theory^, dia¬ 
grammatic approaches [l^ . and the stochastic variational 
method [Ill[l2|- To date, the determination of the energy 
spectrum of systems containing four or more particles re¬ 
mains a challenge, (l^ 

In this paper, we investigate an alternative approach to 
study large systems of fermions. Our symmetry-invariant 
perturbation theory, SPT, offers a perturbation approach 
for a systematic study of correlation including the uni¬ 
tary regime. The perturbation parameter depends on the 
inverse dimensionality of space (6 = I/D), rather than 
the strength of the interaction, so SPT is equally ap¬ 
plicable to weakly or strongly interacting systems. The 
method is essentially analytic [l3j |. with N, the number 
of particles, entering as a parameter, allowing results for 
any N to be obtained from a single calculation[l3. The 
lowest-order result includes correlation and in theory, can 
be systematically improved by going to higher order [l5l|. 
This method couples group theory with conventional di¬ 
mensional perturbation theory (DPT)[I1-[I1 to take ad¬ 
vantage of the high degree of symmetry possible among 
identical particles in higher dimensions. The terms in the 
perturbation expansion are invariant under N\ symmetry 
operations, greatly reducing the number and complex¬ 
ity of the building blocks of these terms. This method 
currently includes full manybody effects that are exact 
through first order. Excited states are obtained from the 
same analytic calculation[l3|. differing only in the num¬ 
ber of quanta in the different normal modes. In principle, 
the full energy spectrum is accessible. 

In past work, DPT has been applied to fermion sys¬ 
tems including single atoms [T^. [T^, quantum dots[^. 
and small molecules (2ll|. The symmetry invariant ap¬ 
proach, SPT, has been applied thus far only to bosons 
with spherical confinement, determining first-order ener¬ 
gies, normal mode frequencies [isl [T3| . the lowest-order 
SPT wave function[2^ and density profilejl^. We ap- 
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plied these results to a BEC for which the density profile 
is a directly observable manifestation of the quantized 
behavior. In a later series of papers, we extended this 
work for bosons to first-order wave functions and density 
profiles [13, [13 ■ We demonstrated that this method effec¬ 
tively rearranges the numerical work for this manybody 
problem into analytic building blocks at each order giv¬ 
ing the exact result order by order in the perturbation 
series [13nil- The complexity of the rearranged problem 
scales with the order of the perturbation series, not with 
the number of oarticlesp^- 

In this extension of SPT to fermions, we tackle the 
challenge of applying the Pauli principle, typically an ex¬ 
pensive numerical task. We describe how this is achieved 
“on paper” for any value of N, thus circumventing heavy 
numerical effort. We apply the method in this initial 
study to large systems of cold fermions in the unitary 
regime for which a number of very accurate calculations 
are available. Results are obtained using analytic build¬ 
ing blocks which have been calculated and stored pre¬ 
viously and which have been extensively checked using 
an independent solution of a model system of harmon¬ 
ically confined, harmonically interacting particles [^. 
Our results compare well to very accurate Monte Carlo 
results H, 1131 including some recent benchmark calcula¬ 
tions using the auxilliary field Monte Carlo method [^. 

We begin the perturbation analysis by defining di¬ 
mensionally scaled quantities: E = k{D)E, and H = 
k{D)H where k{D) is a scale factor which regularizes 
the large-dimension limit . The scaled version of the 
Schrodinger equation becomes 

(^^T+Peff) $ = (1) 

where barred quantities indicate variables in scaled units 
{k{D) = D'^/{ful!ho) for this work. See Ref. [ll]). The 
term T contains the derivative terms of the kinetic energy 
and Raff includes centrifugal, two particle, and confine¬ 
ment potentials [131. 

We assume a totally symmetric, large-dimension con¬ 
figuration at which the effective potential is a minimum. 
The N particles are arranged on a hypersphere, each par¬ 
ticle with a radius, roo from the center of the confining 
potential. Furthermore, the angle cosines between each 
pair of particles take on the same value, 7 ^, i.e. 


ground state, is a zero-point energy contribution satis¬ 
fying the minimum uncertainty principle[30l|- The value 
of 7 oq , which is zero in the mean-field approximation for 
the L = 0 angular momentum states considered here, is, 
in fact, not zero, an indication that beyond-mean-field 
effects are included in the <5^0 limit. 

This highly-symmetric, <5 —>■ 0 structure imparts a 
point group structure to the system which is isomorphic 
to the symmetric group of N identical objects [Slj, Sn, 
allowing a largely analytic solution. The i5 ^ 0 approx¬ 
imation may be systematically improved by u sing it as 
the starting point for a perturbation expansion l^] . The 
Sn symmetry greatly simplifies this task since the inter¬ 
action terms individually have to transform as a scalar 
under the Sn point group. 

The perturbation series has the form: 


E=E^+sY,{s^y 

j=o 

00 


$ = 


E, 




j=o 


(3) 


In practice Ej = 0 \/ j odd. The j = 0 terms are ob¬ 
tained from a harmonic equation, and referred to as the 
energy and wave function at harmonic order. To obtain 
this harmonic correction for small values of 6 , we expand 
about the minimum of the (5—^0 effective potential. 

The harmonic-order Hamiltonian is solved using the 
FG matrix method [s^] to obtain the normal-mode fre¬ 
quencies, uifj,. The number of roots, A^, (A^ = uiff ), of 
the secular equation, N{N l)/2 roots, is potentially 
huge; however, due to the Sj^ symmetry of the problem 
there is a reduction to five distinct roots. 

The FG matrix is invariant under Sn , so it does not 
connect subspaces belonging to different irreducible rep¬ 
resentations (irreps.) of iS'Ar[33- Thus the normal co¬ 
ordinates must transform under irreps. of Sn ■ The 
normal coordinates are linear combinations of the ele¬ 
ments of the internal displacement vectors which trans¬ 
form under reducible representations oi Sn ■ One can 
show that these reduce to two 1-dimensional [N] irreps. 
denoted by 0+,0“, two {N — l)-dimensional [N — 1, 1] 
irreps. denoted by 1+, 1“, and one angular N{N — 3)/2- 
dimensional — 2, 2] irrep. denoted by 2[l3|. 

The energy through harmonic order in S is [l3j | 


lim Ti = Too (1 < * < iV), 

D—^oo ^ 2 ^ 

lim jij = 7 ^ (1 < * < j < iV). 

£)—>-oo 

(This symmetric high-dimensional structure is not unlike 
the localized structure found in a hyperspherical treat¬ 
ment of the confined two-component Fermi gas in the 
iV —>• cxD limit [^.1 In scaled units the <5 —?► 0 (U —J- 00 ) 
approximation for the energy is simply the effective po¬ 
tential minimum, i.e. Fqo = Kff ("^ 00 ) 7 ooi <5 = 0 ). 

In this (5—^0 approximation, the centrifugal-like term 
that appears in I 4 ff, which is nonzero even for the 


E — Eoo + d 


^ (Uf, -k + Vo , 

M={0±1±, J 

2 } 


(4) 


where is the total number of quanta in the normal 
mode with the frequency ^ is a label which runs 
over 0~ , 0+ , 1“ , 1+ , and 2 , regardless of the number 
of p articles in the system (see Refs. [l3j|and Ref. [15] in 
[^1. and Vo is a constant. The multiplicities of the five 
roots are: c?o+ = li do- = 1 , di+ = — 1 , d^- = 

N-1, d2= N{N-3)/2. 
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A character analysis of the normal modes reveals that 
the 2 normal modes are phonon, i.e. compressional 
modes; the modes show single-particle character, and 
the 0^ modes describe center-of-mass and breathing mo¬ 
tions. 

Enforcing the Pauli Principle- To generalize SPT 
from quantum systems of bosons to quantum systems 
of fermions, we must enforce the Pauli principle, thus re¬ 
quiring the A^-body wave function to be totally antisym¬ 
metric. This is enforced by placing certain restrictions on 
the occupancies of the normal modes, i.e. on the values of 
the normal mode quantum numbers, n^, ^ = 0^, 1^, 2 in 
Eq. (|4]) [l^ . The possible assignments can be found by re¬ 
lating the normal mode states |no+,no-,ni+,ni-,n 2 > 
to the states of the confining potential which is a spher¬ 
ically symmetric three dimensional harmonic oscillator 
(Vconf(T'i) = which the restrictions im¬ 

posed by antisymmetry are known. These two series 
of states can be related in the double limit D ^ oo, 
ujho —>• oo where both representations are valid. 

For large D, the normal mode description given by 
Eq. Q is exact. Applying the large uJho limit results in: 


E — M—hiOfiQ {2nQ+ -j- 2nQ- -t- -t- -t- ‘2n2^huJho 

(5) 

Now consider w/jo oo first and then D — ^ cjo. The 
harmonic oscillator levels are exact: 


N 




(2r'i + ^0 + 


hojho 


D 


N 


— N — hlOho + ^ h)ft^ho 


i=l 


( 6 ) 


where is a radial quantum number and U the orbital 
angular momentum quantum number. Equating these 
two expressions which are equal in the double limit, the 
quantum numbers in the two representations can now be 
related to show the restrictions on normal mode states 
imposed by antisymmetry. Because of the clean separa¬ 
tion of radial and angular motions, two conditions result: 


N 

‘2nQ— -j- ‘2n^— — ^ ^ 2i^2 

(7) 

2,TIq+ “h 2/71^+ 2712 ~ ^ ^ li 

These equations determine a set of possible normal mode 
states |71 o+,^o~5^i+:^i-j^ 2 > from the known set of 
permissible L = 0 harmonic oscillator configurations. 

Application: The Unitary Gas- The Schrodinger 
equation for an A^-body system of fermions, N = A^i + A ^2 
with Ni spin up and N 2 spin down fermions, confined by 


E 


70 

60 

50 

40 

30 

20 




o 

° ♦ 





FIG. 1: Ground state energies of the harmonically trapped 
unitary Fermi gas (units Hwho = 1). Our first-order per¬ 
turbation results (filled diamonds) are compared to GFMC 
( + ’8)113, fixed-node DMC (open circles) from Ref. (Hh and 
AFMC results (filled circles, A=6, 8, 14, 20, 30 only) [2^. 


a spherically symmetric potential is 


H'l> = 


N Ni N2 

E ^ E E ^*-5 


'b = E'^ . 


( 8 ) 


where hi and Qij are the single-particle Hamiltonian and 
the two-body interaction potential, respectively. We as¬ 
sume a T = QK condensate with A^i = N 2 confined by 
an isotropic, harmonic trap with frequency uJho- 

To study the unitary regime, we replace the actual 
atom-atom potential by an attractive square well poten¬ 
tial of radius R-. 


Unt {Tij ) 


To , rij ^ R 

0, Ty >R. ’ 


(9) 


For fixed range R, the potential depth Vq is adjusted so 
the s-wave scattering length, Og is infinite. The range 
is selected so R « aho ( aho = \Jh/{moJho)) and can 
be systematically reduced to extrapolate to zero-range 
interaction. We dimensionally continue the square well 
potential so that it is differentiable away from D = 3, 
allowing us to perform the dimensional perturbation 
analysis [T^ fT3| ). Thus, we take the interaction to be 


VintiUj) = Vb(<5) 


1 — tanh 


1 - 3(5 


36R) , 


( 10 ) 

where 14((5) = y+^. The potential depth I 4 is adjusted 
by adjusting the value of b so the scattering length is infi¬ 
nite when 5 = i. This interaction becomes a square well 
of radius R in the physical D = 3 limit. The functional 
form of the potential at I? 7 ^ 3 is not unique. Other forms 
could be chosen with equal success as long as the form 
is differentiable and reduces to a square well potential at 
D = 3. We simply choose a form that allows a gradual 
softening of the square well. 

In Fig. m we plot the ground state energies from 
N = 6 to N = 30 and compare to Green’s function 
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FIG. 2: Ground state energies of the harmonically trapped 
unitary Fermi gas (units hwho = !)• Our first-order perturba¬ 
tion results (filled diamonds) are compared to DMC results 
Ref. [23] (open squares) 

Monte Carlo (GFMC) energies [l^, to fixed-node diffu¬ 
sion Monte Carlo (DMC) energies^ which provide ac¬ 
curate upper bounds to the ground state energy, and to 
recent benchmark auxiliary field Monte Carlo calcula¬ 
tions (AFMC)[2^ which are exact, but subject to finite 
lattice size errors for which corrections have been made. 
The AFMC results are currently the most accurate re¬ 
sults available. Our energies which include full many- 
body effects through first order compare well to these 
AFMC results and except for N = 8 they are slightly 
closer to the AFMC results than the GFMC and DMC 
results shown in Figure 1. Our numbers, as well as the 
GFMC and DMC numbers from Refs. [1, |^, show a dis¬ 
tinct odd/even oscillation, but no obvious shell effects. 
Obtained from analytic building blocks that have been 
calculated and stored previously, our numerical require¬ 
ments take a few seconds on a work station. 

No parameters are used to obtain these results other 
than those used to produce a square well potential with 
infinite scattering length. For the GFMC and DMC stud¬ 


ies shown in Figure 1, no results exist for values of N 
higher than shown. Accurate results for higher N are, 
of course, increasingly difficult to achieve. In Figure 2 
we compare our results at first order above fV = 30 with 
available DMC results from Ref. [2^. Our first-order re¬ 
sults show the expected increase in error as N increases 
suggesting that for larger N, higher order terms may be 
necessary. 

Conclusions - In this paper we have extended the 
symmetry-invariant perturbation method from bosons 
to fermions, applying the Pauli principle “on paper” to 
avoid heavy numerical expense. The method has been 
tested in the unitary regime, which is of particular in¬ 
terest for manybody methods since its infinite scattering 
length and the lack of a natural scale typically require 
intensive numerical simulation for an accurate descrip¬ 
tion. Our analytic results through first order yield ener¬ 
gies that are comparable in accuracy with recent Monte 
Carlo results. As N increases, our error increases sug¬ 
gesting the need for higher order terms. It may also be 
possible to rearrange the perturbation series to minimize 
the importance of higher order terms. 

The theory applied in this paper is applicable to L = 0 
states of spherically confined systems with general at¬ 
tractive or repulsive interparticle interactions and is also 
applicable to both weakly and strongly correlated sys¬ 
tems. The fact that is not zero is an indication that 
beyond-mean-field effects are included in this result even 
in the D —>■ oo limit. This theory is readily generalizable 
to systems with a cylindrical confining potential. 

Higher-order calculations may be required for larger 
N. A detailed program for calculating higher-order DPT 
corrections to A-body systems has been laid out[T^ and 
applied to high order for small-A systems [s^. For large- 
A systems the Sn point-group symmetry greatly simpli¬ 
fies the calculation of these higher-order terms. 

We thank Doerte Blume and Stefano Gandolfi for dis¬ 
cussions and Monte Carlo results. 
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